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ABSTRACT 

A  numerical  scheme  of  the  method  of  characteristics,  the  Hartree 
scheme,  is  developed  to  calculate  the  flow  field  of  the  expansion  of 
a  high  pressure  sphere  into  atmosphere.  The  accuracy  of  this  scheme 
in  solving  one  dimensional  fluid  problems  is  evaluated  by  applying  it 
to  the  blast  wave  problem  for  which  the  exact  solution  can  readily  be 
calculated  to  a  high  degree  of  accuracy.  It  is  shown  that  this  method 
is  very  accurate,  involving  errors  of  less  than  1%.  In  calculating 
the  expanding  sphere,  two  rather  challenging  problems,  namely,  the 
initial  singularity  and  the  formation  of  a  second  shock,  are  suc¬ 
cessfully  solved  through  special  techniques. 

The  propagation  of  the  expanding  main  shock  and  the  propagation 
of  the  contact  surface  which  separates  the  ambient  fluid  from  the 
fluid  initially  confined  in  the  high  pressure  sphere  are  accurately 
determined. 

The  formation  of  an  inward  traveling  shock,  or  the  second  shock 
to  distinguish  it  from  the  main  shock,  from  the  previously  continuous 
flow  field  is  found  to  exist  at  the  tail  of  the  left  traveling 
rarefaction  wave.  Though  the  strength  of  the  second  shock  remains 
rather  weak  at  the  early  stage  of  its  development,  it  grows  rapidly 
to  very  high  strength  just  before  it  implodes  on  the  center  of  the 
sphere.  The  path  of  the  second  shock  in  the  physical  plane  is 
accurately  defined  and  some  interesting  behavior  of  the  second  shock 


is  revealed. 


Based'  on  accurate  detailed  numerical  calculations  by  the 
Hartree  scheme,  it  is  shown  that  "late-stage  equivalence"  exists 
in  the  expansion  of  high  pressure  spheres  into  atmosphere,  provided 
the  initial  total  energy  in  each  of  the  spheres  is  held  constant. 
Late-stage  equivalence  is  assumed  to  exist  if  the  peak  pressure 
distribution  for  different  expanding  spheres  are  the  same  for  long 
times.  Moreover,  if  the  initial  masses  enclosed  in  each  sphere  are 
also  kept  the  same  in  addition  to  constant  initial  energy,  late-stage 
equivalence  exists  riot  only  for  peak  pressure  of  the  main  shock,  but 
also  for  the  positions  of  both  the  main  shock  and  the  second  shock  in 
the  physical  plane. 


NOMENCLATURE 


speed  of  sound 

constant  speed  of  sound  outside  wave  zone 
specific  heat  at  constant  pressure 
specific  heat  at  constant  volume 

parameter  proportional  to  the  total  energy  within  the  wave 
initial  total  energy  released 
enthalpy  per  unit  mass 
constant 

initial  pressure  of  the  high  pressure  sphere 
pressure 

constant  pressure  outside  wave  zone 

entropy  per  unit  mass 

time  variable 

absolute  temperature 

particle  velocity 

shock  wave  velocity 

space  variable 

initial  radius  of  the  high  pressure  sphere 

specific  heat  ratio 

density 

initial  density  of  the  high  pressure  sp.,je 
constant  density  outside  wave  zone 


E 

<j  o 

length  expressing  energy  and  pressure  scaling  (e3  *  — ) 

M 

energy  reduced  radial  distance  (— ' 

clt 

energy  reduced  time  (•*■£—) 


I  INTRODUCTION 


Ever  since  the  discovery  of  explosives  people  have  been  interested 
in  the  behavior  of  explosions.  Today,  considerable  attention  is  still 
focused  on  the  gas-dynamics  of  explosions.  It  is  very  difficult,  if  not 
impossible,  to  represent  the  actual  situation  in  the  initial  stage  for  a 
real  explosion  which  involves  chemical  reactions,  detonation  waves,  com¬ 
bustion  or  nuclear  reactions.  However,  after  the  initial  stage,  the 
situation  is  essentially  the  expansion  of  a  high  pressure  sphere.  Even 
this  problem  of  expanding  gas  has  no  exact  solution.  Simplified  models 
as  well  as  various  numerical  calculations  have  been  used  to  study  the 
problem  in  the  past. 

The  most  well-known  simplified  model  is  the  point  source  solution; 
while  among  the  numerical  calculations  the  finite-difference  method  and 
method  of  characteristics  are  widely  adopted, 

1.  Point  Source 

In  the  so-called  point  source  model  it  is  assumed  that  a  finite 
amount  of  energy  in  an  infinitely  concentrated  form  is  being  released 
suddenly  into  the  atmosphere.  The  solution  of  this  ideal  model  will  be 
referred  to  as  the  point  source  solution  in  later  discussions. 

In  the  point  source  solutions  of  Taylor  ( 1  )  and  Sedov  ( 2  )  the 
properties  of  the  flow  are  assumed  to  be  self-similar;  the  governing 
partial  differential  equations  are  reduced  to  ordinary  differential 
equations. 
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Taylor  Integrated  these  equations  numerically  and  the  results 
vere  presented  in  tabular  form  while  Sedov  successfully  integrated 
these  ordinary  equations  and  obtained  closed  ..ora  solutions.  In  their 
analyses ,  strong  shock  relations  were  used  across  the  outward  propa¬ 
gating  shock  front.  Therefore,  the  solution  is  invalid  for  shock 
pressure  somewhat  below  100  atm.  Besides,,  in  real  explosions  the 
source  of  energy  is  far  from  being  a  point. 

In  order  to  extend  the  range  of  validity  of  the  point  source 
solution,  the  exact  shock  equations  must  be  used  whenever  the  influence 
of  the  atmospheric  counter  pressure  can  no  longer  be  ignored. 

To  extend  to  low  pressures  Goldstine  and  Von  Neumann  ( 3 )  modified 
the  point  source  solution  in  the  range  from  100  atm.  to  1.017  atm.  A 
finite-difference  method  was  used  in  the  calculation  with  initial  data 
calculated  from  point  source  solution.  Exact  shock  relations  were  used 
across  the  shock  front  at  low  pressures.  Numerical  results  were  given 
which  can  be  used  to  compare  with  our  calculation. 

Erode  (4)  also  extended  the  point  source  solution  in  the  range  of 
1,600  atm.  to  1.06  atm.  with  a  finite-difference  method.  However,  Brode 
used  the  artificial  viscosity  method  in  handling  the  shock  discontinuity. 
Since  this  particular  method  always  has  to  spread  the  chock  discontinuity 
across  several  grid  zones,  the  shock  front  cannot  be  clearly  defined. 

In  a  different  approach  Sakural  (5  )  modified  the  point  source 

2 

solution  by  a  series  expansion  in  y,  defined  by  c0,  where  cQ  is  the 

C* 

velocity  of  sound  of  the  atmosphere,  and  D  Is  the  shock  velocity.  The 
seroth  order  solution  is  exactly  Taylor's  solution.  Higher  order 
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approximations  supplement  Taylor’s  solution  such  that  the  exact  shock 
equations  are  satisfied. 

2.  Numerical  Calculations 
A.  Finite  Difference 

One  of  the  earliest  numerical  calculations  by  the  finite-difference 
method  was  by  Unwin  (6),  who  in  1941,  calculated  a  problem  starting 
from  the  sudden  release  of  a  spherically  symmetric  flow  in  a  compressible 
medium,  having  constant  entropy  everywhere,  into  the  atmosphere.  The 
initial  density  distribution  is  so  defined  such  that  no  shock  will  be 
formed.  Results  are  not  applicable  to  the  general  case  where  shocks  are 
involved. 

Roberts  (  7 )  used  the  finite-difference  method  developed  by  P.  D. 

Lax  to  calculate  the  same  problem  treated  by  Unwin.  Lax’s  finite- 
difference  scheme  introduces  a  viscous  effect  which  tends  to  smooth  out 
discontinuities  in  the  flow  variables  when  they  arise.  Since  we  expect 
a  second  shock  which  will  be  formed  in  the  expanding  gas  for  our  problem£ 
therefore,  it  is  reasonable  to  assume  this  method  will  not  yield  satis¬ 
factory  results. 

The  attenuation  of  spherical  shock  at  large  distance  from  the  origin 
is  investigated  by  Whithtm  (8).  As  an  approximation,  the  flow  is 
assumed  to  be  isentropic  in  his  analysis.  He  showed  that  for  a  weak 
explosion,  behind  the  shock  front  an  envelope  of  characteristics  is 
formed  and  hence  «  second  shock  must  appear. 

'  To  handle  hydrodynamic  problems  which  Involve  shock  discontinuities. 
Von  Neumann  and  Richtmyer  (9)  developed  a  numerical  method.  In  applying 


this  method,  an  artificial  viscosity  is  introduced,  and  shocks  are 
s beared  out  so  that  the  mathematical  surfaces  of  discontinuity  are 
replaced  by  layers  of  thickness  comparable  to  several  meshes.  There¬ 
fore,  there  is  a  lack  of  definition  of  the  solution.  Brode  (4)  applied  ' 
this  method  to  calculate  the  spherical  wave  produced  by  sudden  release 
of  high  pressure  gas  into  the  atmosphere  to  very  late  times.  Due  to  the 
nature  of  the  method,  which  cannot  handle  the  initial  singularity  at 
the  outer  boundary,  Brode* s  results  are  inaccurate  even  at  the  early 
stage.  One  of  his  calculations  with  high  initial  gas  density,  the 
overpressure  vs  radius  curve  presented  in  his  paper  is  grossly  in¬ 
accurate. 

B.  Method  of  Characteristics 

Since  the  method  of  characteristics  will  give  correct  treatment  of 
singularities  such  as  centered  rarefactions  and  shock  discontinuities, 
it  generally  leads  to  better  defined  boundaries  and  more  accurate 
results. 

Wecken  (10)  calculated,  by  the  standard  method  of  characteristics, 
the  spherical  wave  produced  by  an  expansion  of  a  high  pressure  gas 
sphere.  Based  on  numerical  grounds  he  predicated  the  existence  of  an 
inward  traveling  second  shock,  which  travels  toward  the  inside  relative 
to  the  moving  gas.  It  will  be  reflected  at  r  *  0  and  then  travel  out¬ 
ward  behind  the  main  shock.  Employing  the  same  method,  Wecken  obtained 
numerical  solution  for  the  explosion  originated  from  detonation  of  a 
spherical  charge,  and  the  results  were  used  to  cospare  experimental 
measurements  obtained  by  Schardin  (11) .  Although  results  showed  overall 


agreement  between  experimental  and  numerical  calculations,  Schardin 
considered  Weckcn’s  calculations  lack  accuracy,  and  recommended  further 
calculations  by  digital  computer. 

The  standard  method  of  characteristics  was  also  used  to  calculate 
the  problem  treated  by  Unwin  {  6)  by  Fox  and  Ralston  (12).  Constant 
entropy  everywhere  was  assumed  in  the  formulation.  Therefore,  their 
solutions  are  restricted  to  problems  where  no  shock  discontinuities  are 
present. 

In  1960,  Zhukov  (13)  using  the  standard  method  of  characteristics 
calculated  a  problem  of  sudden  release  of  a  constant  pressure  sphere 
into  vacuum.  Outside  the  expanding  sphere,  a  spherical  envelope  of 
fluid  exists.  The  primary  Interest  of  this  problem  is  the  interaction 
of  the  main  shock  and  the  outer  fluid  shell.  The  second  shock  was 
formed  as  a  result  of  intersection  of  characteristics  of  one  family 
(II  Characteristics).  After  its  formation  the  second  shock  was  calcu¬ 
lated  for  only  a  short  time.  Details  of  the  numerical  scheme  or  results 
were  not  presented. 

An  alternative  numerical  scheme  of  method  of  characteristics  was 
proposed  by  Hartree  (14).  He  suggested  integrating  equations  of  motion 
along  characteristic  directions  but  Instead  of  using  the  characteristic 
grid,  using  a  given  Interval  in  one  independent  variable.,  It  may  give 
the  results  in  a  more  desirable  form,  particularly  if  the  one  Independent 
variable  happens  to  be  time.  This  method  has  the  advantage  of  yielding 
results  at  a  specific  time,  whereas  if  the  results  are  obtained  on  a 
grid  of  characteristics,  lengthy  interpolation  is  required  to  obtain 


such  results.  Besides,  the  Hartree  scheme  gives  results  along  particle 
paths  thus  making  it  particularly  advantageous  in  handling  problems  with 
contact  discontinuities. 

Hoskin  05)  discussed  Hartree' s  scheme  in  great  detail,  but  no 
specific  example  was  presented  to  demonstrate  its  application. 

Lister  (16)  followed  essentially  Hartree 's  proposal  and  devised  a 
numerical  scheme,  geared  to  digital  computer  computation,  which  can  be 
applied  only  to  isentropic  flow. 

Katshova  and  Chuskin  (17)  developed  a  numerical  scheme,  similar  to 
the  Hartree  scheme,  suitable  for  computer  calculations  of  steady 
axisymmetric  supersonic  flows  of  perfect  gas  with  shock  wave.  Their 
numerical  results  were  in  good  agreement  with  those  by  standard  method 
of  characteristics. 

It  seem3  that,  up  to  the  present,  the  Hartree  scheme  has  not  been 
applied  to  unsteady  flow  problems  Involving  shock  waves  and  entropy 
change. 

One  of  the  diff;  .ulties  in  the  solution  of  the  expansion  of  a  high 
pressure  sphere  is  the  initial  singularity  at  the  boundary.  In  general, 
an  outward  moving  shock,  a  contact  discontinuity  surface,  and  a  centered 
rarefaction  wave  are  formed  instantaneously.  Unwin  (6),  Roberts  (7) 
and  Fox  (12)  avoid  these  difficulties  by  requiring  the  initial  distribu¬ 
tion  of  density  to  be  a  particular  function  of  radius,  and  assuming 
entropy  to  be  constant  everywhere. 

KcFadden  (18)  obtained  a  solution  for  short  times  after  the  release 
of  a  uniform  high  pressure  sphere.  Flow  properties  u,  c,  and  S  (entropy) 
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vere  developed  In  powers  of  y,  which  is  proportional  to  time.  Solutions 
up  to  the  fir3t  order  were  derived  which  satisfy  the  boundary  conditions 
to  the  first  order  without  introducing  a  second  shock.  As  he  suggested , 
we  have  computed  the  example  given  in  his  paper,  and  find  that  our 
results  are  in  good  agreement  with  HcFadden's. 

The  mathematical  nature  of  the  singularity  of  ‘the  disturbance  due 
to  the  detonation  of  a  spherical  charge  initiated  at  the  center  was 
investigated  by  Berry  and  Holt  (19) .  The  equations  of  unsteady  spherical 
motion  are  solved  in  the  neighborhood  of  the  singularity  at  the  origin 
of  the  air  blest  wave  in  the  t-x  plane.  Expansions  are  used  in  series 
of  half-powers  of  r,  the  radial  distance  from  origin,  with  coefficients  0 
depending  on  the  transverse  coordinate. 

Two  singular  characteristics  are  found  to  start  at  the  origin.  One 
is  the  left  traveling  characteristic  in  the  rarefaction  wave  region  and 
the  other  is  the  contact  discontinuity  surface.  Hear  these  lines  the 
expansions  are  invalid.  To  remove  these  difficulties.  Berry  and  Holt  (19) 
used  the  Paincare-Lighthill-Kou  method  by  introducing  a  new  independent 
variable  z  and  expanded  the  dependent  variables  (u,  p,  c)  and  the 
independent  variable  (  Q  and  j  )  in  series  of  power  of  J  with  coef¬ 
ficients  depending  on  z. 

This  is  equivalen.  to  a  perturbation  of  coordinates  Q  .  The 
development  of  the  second  shock  was  investigated  and  asserted  to  be  an 
effect  of  at  least  second  order  of  |  , 

The  results  from  the  series  expansion  were  used  as  initial  data  to 
determine  the  early  development  of  the  growth  of  spherical  blast  from  a 


particular  charge.  Berry  et  ai  (20)  integrated  the  equations  of  motion 
along  characteristics.  Results  of  the  calculation  are  crude.  They  were 
content  with  an  accuracy  adequate  for  plot  only. 

3.  Late -Stage  Equivalence 

It  was  suggested  by  many  authors  (Taylor >  Brode,  Sakurai,  Courant 
and  Friedrichs,  etc.)  that  real  explosions,  though  they  behave  dif¬ 
ferently  in  the  early  stage,  gradually  become  equivalent  if  the  initial 
energies  are  equal.  Among  these  who  suggested  such  a  late-stage 
equivalence,  Sakurai  (5)  made  a  very  strong  statement  that  the  point 
source  model,  although  not  adequate  to  represent  a  real  explosion  in  the 
early  stage,  becomes  more  and  more  accurate  at  the  late  stages,  regard¬ 
less  of  the  kino'  of  explosives  or  the  features  of  the  explosives.  The 
assumption  of  late-stage  equivalence  intuitively  seems  plausible.  How¬ 
ever,  it  has  not  been  verified.  Brode  (4)  calculated  a  few  cases  of  the 
releasing  of  high  pressure  spheres  into  the  atmosphere  and  tried  to 
indicate  that  the  pressure  behind  the  main  shock  front  approaches  the 
corresponding  point  source  solution  in  late  stages.  Unfortunately  due 
to  the  limitation  of  the  numerical  method  used  in  his  calculation, 
results  were  inaccurate  even  at  the  early  stage, 

4.  Outline  of  the  report 

It  is  the  intention  of  this  report  first  to  develop  a  numerical 
scheme  which  is  accurate,  and  can  be  easily  adapted  to  computer  calcu¬ 
lations.  This  scheme  is  then  used  to  calculate  the  flow  field  developed 
from  a  spherically  symmetric  high  pressure  sphere  suddenly  released  in 
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the  atmosphere.  The  existence  of  a  late-stage  equivalence  will  be 
investigated  by  calculating  r  jmparing  expansion  of  many  high 
pressure  spheres  with  different  initial  radii,  pressures,  and  densities. 
Late-stage  equivalence  is  assumed  to  exist  if  the  peak  pressure  distri¬ 
butions  for  different  expanding  spheres  are  the  same  for  long  times. 


II  BASIC  METHOD 


I,  Basic  Differential  Equations 

The  governing  equations  for  an  unsteady  spherically  symmetric 
continuous  flow  of  polytropic  gas  without  friction  and  heat  transfer 
are  the  continuity  equation 

|t  +  2*5t  =  0  (2’1) 

the  momentum  equation 


iU 

TE 


+  u 


ru 

n 


(2-2) 


the  energy  equation 

H  +wtf =o 


and  the  equation  of  state 

5  =  50  *  J5-LOJ&) 


(2-3) 


(2-4) 


where  S0  »  constant 
Substituting  (2-4)  into  (2-3)  yields 


(2-5) 


2.  Method  of  Characteristics  and  Characteristic  Equations 
The  general  theory  of  characteristics  can  be  found  in  many  books, 
such  as  Von  Mises  (21),  Cotirant  and  Friedrichs  (22).  Here  we  will  limit 
our  discussion  to  the  current  problem  which  is  a  quasi-linear  first  order 
system  of  hyperbolic  type  with  two  independent  variables.  The  essense 
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11 


of  the  method  is  to  search  for  a  set  of  characteristic  coordinates.  By 
using  these  characteristic  coordinates,  a  new  system  can  be  obtained 
from  original  partial  differential  equations.  The  differentiation  of 
the  new  system  will  be  considerably  simplified.  For  the  present  case, 
the  differentiations  for  each  characteristic  equation  will  be  along  a 
single  characteristic  coordinate.  Therefore,  numerical  integration 
along  a  characteristic  coordinate  is  particularly  simple. 

A  characteristic  coordinate  is  along  the  tangent  of  a  character¬ 
istic  curve  which  is  defined  as  a  curve  on  which  the  derivatives  of 
the  fluid  properties  are  indeterminate. 

Since  we  will  have  occasion  to  use  the  characteristic  equation  in 
plane  case,  the  derivation  will  be  kept  in  a  more  general  form.  The 
continuity  equation  becomes 


where  >)  ■  1  for  plane  motion,  "0*2  for  cylindrical  motion  and  'J  *  3 
for  spherical  motion;  the  momentum  equation  remains  as: 

4-  1 1  4.  X  —  0  ami's 

a*  e  d*  (2  7) 

and  from  the  equation  of  state  of  a  poly tropic  gas  and  the  definition 
of  sound  velocity 


equation  (2-5)  can  be  written  in  the  form: 


(2-8) 


-  12  - 


For  regions  where  the  flow  properties  u,  p,  and  f  are  continuous, 
three  equations  for  the  total  differentials  du,  dp,  and  dp  may  be 
written.  Combining  these  differential  relations  with  the  governing 
differential  equations  (2-6),  (2-7)  and  (2-8)  we  obtain  the  following 
aet  of  differential  equations  involving  six  partial  derivatives:  , 

*  &X*  f|f»  §£*  an<*  which  are  in  question  for  indetermination. 


Solving  these  equations  for  by  Cramer's  Rule,  we  obtain 

(2-10) 

2>X  P 

where 

dx  dt  0  0  0  0 

o  0  dx  dt  0  0 

O  0  0  0  dx  dt 

D  *> 

p  0  0  0  U.  } 

l 

a  I  T  o  o  o 

o  o  u  1  -uc2  -c2 

a  (  uJt-Jx)  [(ui-c)Jt  -dx]  [(u-c)dt-dx] 


—  * 
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and 


du 

dt 

0 

0 

0 

0 

0 

dx 

dt 

0 

0 

d? 

0 

0 

0 

dx 

dt 

0H)f 

0 

0 

0 

li 

1 

0 

1 

\ 

J  0 

0 

0 

0 

0 

u. 

1 

-ucz 

-c? 

•*  (ad' t-dx) [(>)-i)-^-(dt)  —  ( udt—dx)du  +  ^-^-dt J 


In  order  that  ^  to  be  indeterminate  D  must  be  0. 

produces  the  three  physical  characteristics: 

The  vanishing  of  D 

I  characteristic 

at  =  utC 

(2-11) 

II  characteristic 

&  =  tt-c 

(2-12) 

III  characteristic 

at  =  w 

(2-13) 

Notice  that  the  III  characteristic  is  a  path  line. 

Along  these  physical  characteristics,  the  derivative  is 

2a 

>U 

Indeterminate  and  therefore  may  be  discontinuous.  To  insure  that 
is  indeterminate  but  not  infinite,  the  numerator,  of  equation  (2-10) 
must  also  vanish  along  the  characteristics.  Along  the  III  character¬ 
istics,  (udt  -  dx)  *  0,  N  *  0  identically;  the  vanishing  of  N  along  the 
other  two  characteristics  yields  the  following  state  characteristics, 

(or  compatibility  equations) 

Wu.  =  ?  ;  Co-O-Sf-Jt  <2-«> 
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where  the  upper  sign  refers  to  I  characteristic  and  the  lower  sign 
refers  to  II  characteristic. 

Wien  equations  (2-3)  are  solved  for  the  other  derivatives,  we  have: 


=  -Jx)  [-($'0  dxdt~(  udu+  ^r)dx+  (u2-C*)dadtl 


(2-15) 


^  =  ^(uWt-clx)^|-(0-0  Qjf-dt  ~df>J(udt-dx)+  fc*dudtj  (2-16) 

£~dx+ud jj (uJt-Jx)~C*d udxj  (2-17) 

The  vanishing  of  the  numerators  of  ^  and  yields 

results  which  are  identical  to  those  for  -|£  . 

Vnen  equations  (2-9)  are  solved  for  and  we  have 


^  ^(u.dt-dx)[(-(-0~t)  Zjfdt  -  cl  0  fudt-dx) 

+  p Judtl  i-  C*d^(dt)*- 

^  udt-dx)£((*-0££dx  +  ud^(udt-dx) 

-  ^Jxduj  +Jf>dxdt  -  UC *c/^(cltfj 


(2-18) 


(2-19) 


The  vanishing  of  the  numerators  of  and  yields,  in  addition 

dX  o V 


to  equation  (2-14),  the  III  state  characteristic 


=  ~CrJf> 


(2-20) 


It  is  convenient  to  eliminate  the  density  f  from  the  characteristic 

equations  by  the  relation  C*~Y"r-'  When  this  substitution  is  made,  and 

K 


3 


I 


1 

I 

I 

I 

1 

I 

I 

I 

I 

1 

I 

I 

1 

1 

I 

I 

I 

1 
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the  III  state  characteristic  equation  integrated,  the  following 
characteristic  equations,  in  terms  of  p,  u,  and  c  are  obtained. 
Stite  Characteristics 


-  -  jf  (Jf)t  -  ('>-') 

(2-21) 

(2-22) 

(fi = ©f 

(2-23) 

The  constants  c*  and  p*  must  be  evaluated  for  each  III  character¬ 
istic. 

In  the  numerical  procedure  of  the  method  of  characteristics,  all 
governing  differential  equations  are  put  into  finite-difference  form. 
Thus,  equations  (2-11)  to  (2-13),  (2-21)  and  (2-22)  become 


along  I  and  II 

AX 

At 

=  u  ±  C 

(2-24) 

(2-25) 

along  III 

AX 

aT 

=  U 

(2-26) 

where  a  represents  the  difference  between  two  adjacent  points  along  a 
characteristic  and  the  barred  values  represent  the  average  between  the 
two  points.  Notice  that  referring  to  the  pivotal  point  (mesh  point) 
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Xj  j,  AU-  IX-  —  U.-_|  is  the  backward  difference  of  u;  although  it  is 
sometimes  considered  as  the  central  difference  at  "half-way"  point 

*■  *  • 

3.  Normal  Shock  Equations 

The  equations  expressing  the  conservations  of  mass,  momentum  and 
energy  across  a  moving  shock  are: 


e,n  =  e*  [  v  t 

(  U  ru2 )] 

(2-27) 

P,  +  P,  u2  =  f. 

*  ?3[U  t  (u,-u2)] 

(2-28) 

H,  +  T  = 

.  [Vt(u,-u,)l2 

2 

(2-29) 

where  U  is  the  shock  velocity  relative  to  the  medium  into  which  it  moves. 
Subscripts  1  and  2  refer  to  the  states  ahead  of  and  behind  the  shock 
front  respectively.  The  upper  sign  applied  for  shocks  travels  to  the 
right  and  the  lower  sign  for  left  traveling  shocks,  but  U  is  always 
taken  as  positive.  To  make  these  equations  more  suitable  for  later  use, 
the  following  algebraic  manipulations  are  in  order.  For  polytropic 
gases  we  have 

H  =  cfr  =  cr  Mf)  (f)~  FTjff) 

Since 

*(■£) 

c* 

therefore  H  -  '  (2-30) 


Using  relation  (2-30) ,  equation  (2-29)  becomes 


c/-c/=  -  "T-r  u,-u2)[(uru2)  t  2  u] 


From  equation  (2-27) ,  we  get 


P  .  _ 

V  t  (  u.-a,) 


Substituting  (2-32)  in  (2-28)  we  obtain 


ft  -  ft  =  i  f.UCu-Uj) 


Equation  (2-33)  can  be  written  as 


=  ±  e,V(u--ut) 
r  v  v 


(2-31 


(2-32; 


(2-33] 


(2-34; 


Eliminating  P2  from  (2-34)  with  equation  (2-32)  and  rearranging,  we  get 


c; - 


u  cf 


=  ±  Y  V( 


(2-35; 


2 

Elimination  of  c,  from  (2-31)  with  the  help  of  equation  (2-35) ,  we  have 


U2  -  U,  - 


C  /  1/  £l) 

mU 


(2-36] 


Substituting  u«  from  equation  (2-36)  into  equation  (2-33),  writing 


1 


f»=  £n<*  t^ien  solving  ^or  ^  yields 

v  r  c'JWUn~~0 + 1 


(2-37) 


Substituting  Uj  from  equation  (2-36)  into  equation  (2-31)  and  solving 
for  yields 


=  Ci/J  I  + 


(2-38) 


Equations  (2-36),  (2-37)  and  (2-38)  are  the  ones  actually  used  in  the 

calculations.  For  very  strong  shocks  where  -ik»  j  ,  equation  (2-37) 

Pi 

can  sometimes  be  approximated  by 


v  * 


(2-39) 


It  is  evident  from  equation  (2-39)  that  -£-»/.  Consequently 
equations  (2-36)  and  (2-38)  can  be  approximated  by  the  following  simple 
forms: 


u2  =  u.  ±  J+J 


(2-40) 


'lllilzll  TJ 
S?  ~  ?+i  U 


(2-41) 


Strong  shock  equations,  (2-39),  (2-40)  and  (2-41)  are  used  in  the 
calculations  of  the  blast  waves  where  these  relations  are  consistent 
with  the  similarity  solution  of  the  blast  wave. 


Ill  SOLUTION  OF  THE  PRESENT  PROBLEM 


1.  Choice  of  the  Basic  Numerical  Technique 

There  are  two  major  basic  numerical  techniques  in  the  application 
of  the  method  of  characteristics.  One  is  the  standard  technique  of 
integration  along  the  two  main  characteristics,  while  the  other  follows 
the  Hartree  scheme,  which  involves  fixed  time  intervals.  These  techniques 
were  discussed  and  their  applications  were  demonstrated  in  (23) ,  (24)  and 
(25) ,  where  it  has  been  shown  that  both  numerical  schemes  produce  very 
accurate  results  when  applied  to  the  blast  wave  problem,  provided  the 
region  with  extremely  high  speed  of  sound  is  excluded.  It  is  reasonable 
to  assume  that  either  of  the  two  methods  will  yield  accurate  results  when 
applied  to  the  expansion  of  high  pressure  sphere. 

In  applying  the  standard  technique,  the  points  of  intersection  of 
two  families  of  characteristics  do  not  occur  at  equal  intervals  in  either 
the  space  or  time  coordinates.  When  the  spatial  distribution  of  the 
flow  variables  at  a  later  instant  is  required,  we  have  to  perform  a 
two-dimensional  interpolation,  which  is  not  always  desirable. 

The  constant  time  technique  overcomes  the  difficulty  of  uneven 
distribution  of  grid  points  by  choosing  grid  points  on  predetermined 
constant  time  lines.  This  has  the  advantage  that  the  required  inter¬ 
polation  is  always  one  dimensional. 

The  Hartree  scheme  gives  results  along  particle  paths  in  successive 
steps,  therefore  it  is  advantageous  in  dealing  uith  present  problem 
where  an  interface  (contact  surface)  is  involved.  Difficulty  may  arise 
in  applying  the  standard  technique  to  solve  problems,  where  shock 
formation  is  expected.  Suppose  the  solution  along  characteristic  OA 
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and  to  the  left  of  it  (see  Fig.  1)  has  been  obtained  without  the  prior 
knowledge  of  the  formation  of  a  shock  at  B  which  is  the  intersection  of 
two  characteristics *of  the  same  family.  Let  BL  be  the  subsequent  shock 
path;  then  the  solution  in  the  region  ABL  is  wrong,  since  no  account  was 
taken  of  the  shock  to  establish  it.  To  obtain  the  solution  in  this 
region,  the  data  required  such  as  solutions  at  points  1,  2,  3  and  4  are 
generally  destroyed  in  order  to  save  storage  space.  This  difficulty  can 
be  overcome  either  by  keeping  the  points,  at  which  the  shock  is  most  likely 
to  occur,  at  greater  time  values  than  their  neighboring  points  or  by  retain¬ 
ing  large  amounts  of  data,  which  might  be  useful  for  corrections  due  to  the 
occurrence  of  the  shock,  in  the  memory.  The  former  method  causes  compli¬ 
cations  in  programming  logic  and  the  latter  needs  extra  memory  spaces. 

They  are  both  undesirable  in  computer  calculations. 

The  shortcomings  of  the  Hartree  scheme  are  generally  believed  to  be 
that  it  is  less  accurate  and  more  time  consuming.  Nevertheless,  the 
results  of  the  calculation  for  the  plane  blast  wave  by  the  standard 
technique  and  Hartree  technique  showed  that  the  same  accuracy  was  reached 
by  both  methods  within  comparable  computer  time.  The  main  reason  for  this 
favorable  result  is  because  in  applying  the  constant  time  method,  the 
initial  data  points  are  on  the  same  time  line,  therefore  quadratic  inter¬ 
polation  can  easily  be  used  in  determining  the  flow  properties  at  points 
other  than  the  initial  data  points.  As  quadratic  interpolations  are 
generally  more  accurate  than  linear  interpolations,  this  in  turn  improves 
the  accuracy  of  the  new  point  to  be  calculated.  Consequently,  in  applying 
the  Hartree  technique,  we  are  able  to  start  with  fewer  initial  data  points 
in  the  beginning  of  calculation  and  still  attain  the  same  degree  of  accuracy. 
As  a  result,  this  compensates  to  the  same  extent  the  time  required  in  carrying 
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out  extra  interpolations.  It  may  also  be  attributed  to  the  more  even 
distribution  of  data  points  throughout  the  region  of  calculation  when  the 
Hartree  scheme  is  used.  By  now  it  is  evident  that  the  Hartree  technique 
is  preferable  for  application  to  the  present  problem. 

Although  the  basic  numerical  technique  for  the  Hartree  scheme  of 
method  of  characteristics  has  been  developed  in  reference  (25),  nevertheless, 
due  to  the  complexity  of  the  present  problem,  quite  a  few  difficult  points 
still  need  to  be  resolved.  In  this  section,  the  method  of  handling  the 
initial  singularity,  methods  of  initiation  of  the  second  shock,  and  some 
other  problems  pertaining  to  the  solution  of  the  present  problem  will  be 
discussed. 

2.  Starting  Singularity 

Immediately  after  the  release  of  the  high  pressure  sphere,  an  outward 
traveling  shock,  which  will  be  referred  to  later  as  the  main  shock  to  dis¬ 
tinguish  it  from  the  second  shock  formed  at  a  later  time,  propagates  into 
the  stagnant  air;  in  addition,  a  contact  surface  which  separates  the  gas 
and  the  air,  and  a  centered  rarefaction  wave,  are  initiated  from  the 
boundary.  (See  Fig.  2)  The  existence  of  a  second  shock  at  the  tail  of 
the  rarefaction  has  been  discussed  by  many  authors,  i.e.,  Wecken  (10), 

Whit’nam  (8),  Berry  and  Holt  (19),  etc.  It  is  also  known  that  the  strength 
of  such  a  shock  is  zero  at  the  beginning  and  remains  very  weak  in  the  early 
stage  of  the  expansion.  (See  initiation  of  second  shock.)  Therefore,  to 
avoid  further  complication  cf  the  already  very  intricate  starting  singularity, 
no  second  shock  was  introduced  in  the  initial  singularity.  This  introduces 
practically  no  error  at  all,  since  the  error  in  not  considering  the 
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second  shock  is  of  third  order,  which  is  practically  zero  at  this  stage. 

At  t  *  0+,  the  time  increment  At  approaches  to  zero  and  dif¬ 
ferences  in  characteristic  equations  between  the  plane  and  spherical 
cases  vanish. 

For  a  given  high  pressure  sphere  surrounded  by  atmospheric  con¬ 
ditions  all  the  flow  properties  ere  specified  for  points  1  and  2.  What 
has  to  be  established  are  the  properties  at  points  3  and  4.  Point  3, 
situated  between  the  tall  of  the  rarefaction  wave  and  the  contact  sur¬ 
face  represents  the  expanded  gas  on  the  left  of  the  contact  surface, 
while  point  4  corresponds  to  the  air  engulfed  by  the  outward  propa¬ 
gating  shock  on  the  other  side  of  the  contact  surface.  Notice  that 
across  the  contact  surface,  the  pressure  and  velocity  are  continuous, 
i.e.,  Pj  ■  p^  and  «  u^.  However,  Che  speed  of  sound  c^  nay  be 
different  from  c^.  This  depends  on  the  initial  conditions  of  the  high 
pressure  gas  and  that  of  the  atmospheric  air. 

Flow  properties  at  points  3  end  4  can  be  determined  analytically 
fey  the  shock  relationships  'jetween  points  1  and  4,  the  isentropic 
expansion  relation  between  points  2  and  3,  and  the  boundary  conditions 
across  the  contact  surface,  u^  *  u^  and  »  p^. 

For  the  centered  rarefaction  wave  the  flow  is  isentropic,  and 
points  2  and  3  are  on  a  I  characteristic  line,  When  the  state 
characteristic  equation  (2-14)  is  simplified  for  the  plane  case  (  'i  »  1), 
we  have 

(<K  =  - 


(3-D 
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For  isentropie  flow  of  polytropic  gas  with  t  «  Y  the  following 
relationships  hold 


(3-2) 


=»  constant  * 


(3-3) 


With  equations  (3-2)  and  (3-3),  we  can  write: 

c!=  / f/r2frj'(  (3-4) 

Putting  equation  (3-4)  in  differential  form,  wa  obtain 

2  Cdc  ~  K,  £  (r3~l)  f  2  df  (3-5) 

2 

Dividing  (3-5)  by  c  ,  we  get 

2 =  (<i -O^f  0-6) 

Differentiating  (3-3)  and  rearranging  yields 

dp  =  =  C2d^>  (3.7) 


Eliminating  dp  from  equations  (3-7)  and  (3-5)  we  have  the  expression: 


ZC  =  c r*-0  j£- 


(3-8) 
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Again  by  elimination  of  dp  from  equation  (3-1)  and  (3-8),  we  finally 
get 


(du)i  =  -  jpj  (dc\ 


(3-9) 


Therefore 


u3  '  u2  “  -7^jCc3-C2).  With  u2  =  0 


W?  =  f,-  f  (  \ 3  ~  ) 


This  equation  can  be  written  as 


^  =  2  _  _£l\ 

c,  ^Tl'  cj 


(3-10) 


Introducing  the  isentropic  relation  between  points  3  and  2 


equation  (3-10)  becomes 


(3-11) 


For  the  right -vard  traveling  shock  with  u-  “0  and  X  =  X,  ,  the  cor¬ 
responding  shock  equations  (2-36)  and  (2-37)  are 


U 


4- 


_2_(  V* 

r,+  i  l 


A' 


(3-12) 
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and 


(3  -13) 


Eliminating  from  equations  (3-12)  and  (3-13),  we  obtain 


(3-14) 


Since  u^  *  u3  and  p^  -  p3,  equation  (3-11)  is  equivalent  to 


(3-15) 


Equating  the  expressions  for  u^  in  equations  (3-14)  and  (3-15)  and 

Y 

solving  for  -jr  ,  we  get: 


(3-16) 


Here  in  equation  (3-16),  p^  is  the  only  unknown.  Once  the  initial 
conditions  of  the  gas  as  well  as  that  of  the  atmospheric  air  are 
known,  can  be  evaluated. 

The  Newt&n-Raphson  iterative  process  was  used  in  the  calculation 


of  and  no  difficulties  have  been  encountered.  After  p^  is 
determined,  the  remaining  unknowns  can  be  obtained  in  the  following 
fashion:  the  shock  velocity  can  be  calculated  from  equation  (3-13), 


JL 

then  and  can  be  evaluated  from  equations  (3-12)  and  (2-38) , 
respectively.  Since  =  p^  and  s  u^,  the  only  unknown  left  is 
c^  which  can  easily  be  obtained  from  equation  (3-10).  Now  everything 

J. 

has  been  determined  for  x  «  x  and  t  **  0  . 

o 

Again,  referring  to  Fig.  2,  at  time  t  “  tj  “  At  (a  small  time 
interval)  the  physical  locations  of  points  8,  N,  5  and  7  as  well  as 
the  flow  properties  between  8  and  7  are  to  be  determined. 

For  a  very  small  time  interval,  At,  the  flow  properties  in  the 
regions  N05  and  607  vary  very  little.  Therefore,  in  changing  the 
differential  equations  into  corresponding  finite  difference  equations 
and  interpolating  properties  between  points  can  be  made  without 
incurring  much  inaccuracy.  However,  this  is  not  the  case  for  region 
80N,  the  centered  rarefaction  wave  region,  where  flow  properties, 
even  though  continuous,  vary  very  rapidly  with  respect  to  both  x  and  t, 
no  matter  how  small  At  is.  To  overcome  this  difficulty,  the  centered 
rarefaction  wave  region  was  divided  into  many  small  segments  separated 
by  the  II  characteristics  issuing  from  point  0. 

At  point  0,  the  pressure  difference  between  neighboring  II 
characteristics  are  chosen  to  be  equal.  It  is  known  that  flow  proper¬ 
ties  change  gradually  along  each  II  characteristic.  Therefore  the 
variation  of  flow  properties  within  each  segment  will  be  small.  Hence 
it  is  permissible  to  obtain  the  solutions  for  points  from  8  to  N 
successively  with  finite  difference  equations  along  characteristic 
directions. 


*  Replacing  subscript  2  with  4, 
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Since  the  head  of  the  rarefaction  wave  propagates  with  a 

velocity  u  -  c,  where  u  ■  0  and  c  *  C2  “  constant  for  the  undisturbed 

gas,  therefore,  08  traces  a  straight  line  with  slope  -  -C2. 

Knowing  the  location  of  point  8  and  all  the  flow  properties  along 

28,  the  solution  for  point  9  can  be  established. 

Let  9*  be  a  point  on  09  near  point  0,  then  the  pressure  at  point 

9‘  fV  =  +>  JLlJL,  where  K  is  the  total  number  of  segments  chosen 

M 

in  the  rarefaction  wave  region.  The  particle  velocity  u^,  and  speed 
of  sound  Cgt  can  be  calculated  from  equations  (2-23)  and  (3-9), 


where  u^  *  0 

Knowing  the  location  of  point  8,  all  the  flow  properties  along 
28,  and  the  initial  slope  of  the  II  characteristic  line  at  9',  the 
solution  for  point  9  can  be  established  with  the  help  of  the  six  finite- 
difference  equations  along  the  three  characteristic  directions.  This 
is  accomplished  by  iterative  process  similar  to  the  one  used  for  the 
interior  points,  (see  Ref.  (25))  with  the  only  exception  that  the 
initial  data  are  not  located  at  a  constant  time  line. 

Essentially  the 'same  procedure  can  be  used  to  obtain  solution  for 
points  10  to  N  successively.  Here  the  distribution  of  flow  properties 
along  each  of  the  II  characteristics  is  assuiued  to  be  linear. 
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The  solutions  for  points  5,  6  and  7  are  interrelated.  They  have 
to  be  solved  simultaneously  through  a  rather  complicated  system  of 
algebraic  finite-difference  equations.  They  are  (refer  to  Fig.  2) 
finite-difference  equations  derived  from  physical  characteristic 
equations , 


Along  I-  ~ 

8 


u4 » “s  .  (*  +  <* 

Z  2. 


(3-17) 


Along  IIIg 


U3  +■  Us 


(3-18) 


*7-*c  =  Ut  +  u7  ^  Cc  +  C7 


Along  Ia  2 


(3-19) 


^<•“^5  Ug  Cg+'Ctf 

.Along  II.  "  —2 - 2 


(3-20) 


Along  III.  dt 


(This  equation  is  identical 
with  equation  (3-18)  because 

*6"*5»  u4”u3  and 


And  finite-difference  equations  derived  from  state  characteristic 
equations : 
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(3-22) 


Along  IIIg 


Along  I,  u, -uc  =  (3-M) 


Along  IIa  «*-«,-  1f^RyCfc-W  +  ^^(ot-t,)  (3-24) 


Along  IIIa 


(3-25) 


The  shock  relations 


(3-26) 


*7 


(3-27) 


c7 


(3-28) 


By  the  definition  of  velocity  of  the  shock,  we  can  write 


(£x\  _  If  =  Hz  *  ^  ~ 

'^'shock  ^ 


<3»29) 
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Assume  that  At  is  small  enough  such  that  within  the  region  under 
consideration,  the  tail  of  the  rarefaction  wave,  the  contact  surface 
and  the  shock  front  can  be  approximated  by  straight  lines.  And  further¬ 
more  we  assume  that  the  flow  properties  are  linearly  distributed  along 
these  lines.  Then  we  may  write 

At  ~  tA  (3-30) 


__  Xc-*0 

At  “  tc  (3-31) 


% o  _  ~  * o 

At  ”  t3  (3-32) 


and  nine  linear  interpolation  equations  for  flow  properties  at  points 
A,  B  and  C,  Since  solutions  at  points  3  and  4  are  known  from  the 
previous  calculations,  suppose  At  has  been  preselected  (to  be  dis¬ 
cussed  later) ,  then  the  total  number  of  variables  involved  in  the  . 
above  equations  are  28,  namely:  %s  ur  ce  t>„  y  ^ 

>  ■*./  >5,  KXt)}  Ybt  *7  , 

W7  ,  c7,  U7,  UA>  CA>  P/»,  t5j  UBt  CB)  tC)  uCf  Cc  and  f>c  . 

Because  5  and  6  are  points  on  two  sides  of  the  contact  surface,  ur  »*  u... 

5  o’ 

P5  "  Pg  x5  "  xg*  tn*s  reduces  the  total  number  of  unknowns  to  25. 
Equations  (3-17)  to  (3-32)  comprise  16  equations.  Adding  nine  inter¬ 
polation  equations,  we  have  exactly  25  equations  to  solve  for  the  25 
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unknowns.  The  solution  o£  the  above  system  of  equations  was 
accomplished  by  iterative  processes. 

The  accuracy  of  the  numerical  scheme  was  checked  by  applying  it 
to  the  numerical  exanq>le  given  by  McFadden  (18)  and  the  results  were 
compared. 

McFadden  took  the  case  of  the  expansion  of  a  high  pressure  gas 
sphere,  initially  at  rest,  into  surrounding  air.  Both  the  high 
pressure  gas  and  the  ambient  air  are  assumed  to  be  diatomic  gases 
with  t  «  1.4.  The  initial  pressure  ratio  of  the  gas  to  that  of  the 
surrounding  air  is  12.8173  and  the  density  ratio  is  3.9560.  Com- 
parable  results  at  a  time  when  the  head  of  the  rarefaction  wave 
traveled  a  distance  of  5%  of  the  initial  radius  of  the  gas  sphere 
are  as  follows. 


McFadden' s  series 

Results  by  method  of 

solution 

characteristics 

V 

1.055  t 

1.055 

xi 

1.035  t 

1.035 

xt 

0.992  t 

0.992 

* 

4.33  ** 

4.333 

li 

fl 

4.29  ** 

4.298 

* 

4.24  ** 

4.236 

*  x  denotes  the  physical  position  and  the  subscripts  s,  i,  and  t 
are  used  for  the  shock  front,  the  contact  surface  and  the  tall 
of  the  rarefaction  wave,  respectively, 
t  values  calculated  according  to  McFadden' s  solutions. 


**  values  presented  by  McFadden. 

Considering  the  fact  that  these  results  by  method  of  character¬ 
istics  were  obtained  through  a  single  time  step,  it  is  gratifying  to 
see  such  a  good  agreement  with  McFadden* a  series  solution.  In  the 
actual  calculation  of  the  high  pressure  sphere,  smaller  time  steps 
(approximately  half)  are  employed.  Therefore,  the  claims  for  accuracy 
of  this  numerical  procedure  are  warranted. 

3.  Second  Shock 

After  the  release  of  a  high  pressure  sphere,  in  addition  to  the 
main  shock  propagating  into  the  atmosphere,  there  exists  another  shock 
wave,  which  is  usually  referred  to  as  the  second  shock,  formed  at  the 
tail  of  the  rarefaction  wave  traveling  into  the  expanding  gas.  The 
existence  of  the  second  shock  was  first  predicted  by  Wecken  (10)  on  the 
basis  of  extensive  numerical  calculations  of  a  spherical  blast  from 
certain  explosions.  Berry  and  Holt  (19)  analyzed  the  initial  disturbance 
of  the  spherical  blast  from  certain  explosions  by  series  expansions 
and  showed  that  the  singular  characteristic  initially  in  the  direction 
of  the  tail  of  the  rarefaction  wave  developed  into  a  shoek  wave.  The 
same  phenomenon  is  found  to  occur  for  the  expansion  of  a  uniform  high 
pressure  sphere  as  a  result  of  the  formation  of  a  compression  wave 
behind  the  tail  of  the  rarefaction  wave.  If  a  very  small  left  traveling 
discontinuity  is  introduced,  it  grows  steadily  into  a  finite  shock. 

The  second  shock,  though  growing  with  time,  remains  very  weak 
until  the  time  when  the  head  of  the  rarefaction  wave  reaches  the  center 
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of  the  sphere,  then  it  starts  to  grow  in  strength  rapidly.  The  strength 
of  a  shock  is  defined  as  the  ratio  of  the  difference  in  pressure  across 
the  shock  to  the  pressure  in  front  of  the  shock.  Since  there  is  no 
second  6hock  at  t  »  0,  and  the  initial  rate  of  growth  is  slow,  there¬ 
fore  exactly  where  or  how  large  an  initial  shock  discontinuity  to  be 
introduced  in  the  beginning  is  of  little  influence  in  the  general 
results,  as  long  as  it  is  inserted  before  the  head  of  the  rarefaction 
wave  reaches  the  center  and  the  magnitude  of  the  discontinuity  is  very 
small.  This  fact  can  be  observed  in  Fig.  3  and  Fig.  4,  where  curves 
showing  position  and  the  strength  of  second  shock  vs  time  are  plotted 
for  shocks  initiated  at  different  times  and  with  different  initial 
discontinuities.  These  curves  also  can  serve  the  purpose  of  showing 

the  general  behavior  of  the  second  shock.  It  is  called  a  left  traveling 

shock  because  it  propagates  toward  the  left  with  respect  to  the  gas 
particles  in  front  of  it.  Depending  on  the  values  of  the  particle 
velocity  u  and  the  shock  velocity  U  ,  the  absolute  velocity  of  the 
shock  front  may  move  toward  either  left  or  right.  However,  as  the 
shock  grows  in  strength,  it  will  eventually 'acquire  an  absolute 
velocity  toward  the  center  of  the  sphere. 

The  general  procedure  in  the  initiation  of  the  second  shock  is 
described  as  follows: 

Suppose  a  second  shock  is  to  be  initiated  at  the  time  t  ■  t2  (see 
Fig.  5 ) ,  then  the  calculation  in  the  rarefaction  wave  region  will  end 
at  point  8*  which  is  on  the  II  characteristic  immediately  adjacent  to 

the  left  of  the  tail  of  the  rarefaction  wave. 
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The  chock  at  point  7  is  initiated  by  inserting  a  s'uock 
discontinuity  of  negligible  strength  at  point  2  on  the  initial  time 
line.  Let  .Q~,.fAr g  «  then  u^,  and  can  be  calculated  from 
the  shock  relations 


Now  everything  is  known  for  points  8r  1,  2,  3,  4  and  5.  Write  the 
finite-difference  equation  for  physical  and  state  characteristic 
equations . 

Along  I^-characteristic 

X^-X/i  Va  +  uc,  4.  CA  +  C6 

~Vt„  2  2  (3-33) 


Ufc-tU 


C(,+  CA 

\ThFh) 


LK.fAyh0^^A) 


(3-34) 


Along  Illj-eharacteristic 

Xfe  ~  Xx>  _  Uj>  +  Uu 

~  "t*  2 

A  = 

fp  '■C,/ 


(3-35) 


(3-36) 
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Along  IIj -characteristic 

*6"  *8  _  Kgt-tffc  £3—37) 

“  At  "  ~Z  2 


u6-u8  =-  *VCb  .  (p,-K)  4-  At 

Along  II, -characteristic 


(3-38) 


*6-*:  -  U7*-Ut  _  C7+Cc 

/Tt - 2 —  — 2~ 


(3-39) 


<7 +  <C  h  ^7^)^  4*0;)  At 

“7  '  %-(h*T3(h  r‘} 


(3-40) 


Across  the  shock  we  can  write: 


(3-41) 


(3-42) 


C 


7 


2(r*-0 

(«*'/ 


(3-43) 


From  the  definition  of  shock  velocity 

-f- -  U7 

TE  2 


(3-44) 
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Assuming  the  11^  characteristic  18  is  a  straight  line,  then  two  more 
equations  can  be  written: 

_At  _  t,-tA 

V*i  ~  *b-Xa  <3~*5> 


A  t  K  ~  ~tj, 
'Xg'Xt  P 


(3-46) 


By  assuming  the  flow  properties  are  linear  along  18  and  12  yields 

nine  equations  for  linear  interpolation  of  properties  at  A,  D  and  B. 

Since  the  flow  properties  varying  rapidly  between  points  3,  4  and  5, 

three  quadratic  interpolation  formulas  were  used  in  determining  u  ,  c  , 

C  C 

and  p  .  Adding  the  12  equations  from  interpolating  the  flow  properties 

V 

to  the  14  equations  (3-33)  to  (3-46),  makes  a  total  of  26  equations. 

There  are  exactly  26  unknowns,  namely:  X5,  Ug,  c^,  p^,  U7,  u7,  c 7,  p7, 

JtA*  fcA»  UA>  CA»  PA»  XD»  fcD»  UD»  CD‘  PD*  *B»  UB»  CB»  PB»  XC»  UC»  CC’  and 
p  .  The  approximate  solution  of  this  system  of  equations  can  be  arrived 

V* 

at  through  an  iterative  process  without  difficulty, 

4,  Contact  Surface 

One  of  the  advantages  of  applying  the  Hartree  scheme  in  the  calcu¬ 
lation  of  the  present  problem  can  be  appreciated  from  the  easiness  in 
treating  the  contact  surface.  Since  the  contact  surface  it?  formed  by 
particle  paths,  essentially  the  same  numerical  scheme  as  for  the  interior 
points  can  be  used  with  the  only  difference  being  that  the  isentropic 
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relations  along  particle  paths  yield  two  equations  rather  than  one 
equation  in  the  case  of  an  interior  point.  This  is  because  across  a 
contact  surface  pressure  and  particle  velocity  are  continuous  but  the 
speeds  of  sound  usually  are  different.  Finite-difference  equations 
along  characteristic  directions  are  similar  to  those  of  the  interior 
points  which  are  shown  in  reference  (25).  In  the  solution  of  these 
equations  almost  the  same  Iterative  procedure  for  the  interior  points 
can  be  employed. 

5,  Choice  of  Time  Step 

The  choice  of  the  magnitude  of  the  first  time  step  At^  is 
arbitrary  except  for  reasons  of  accuracy  and  computational  effort. 

This  time  step  is  usually  kept  in  a  range  that  in  such  a  time  interval 
the  head  of  the  rarefaction  travels  a  distance  approximately  equal  to 
3Z  of  the  initial  radius  of  the  high  pressure  sphere. 

a  V 

In  the  subsequent  calculations,  the  Courant  condition  At  <  -jr-, 
which  was  discussed  in  reference  (25),  is  observed.  However,  before  the 
head  of  the  rarefaction  reaches  the  center,  the  choice  of  time  steps 
are  further  limited  to  be  not  greater  than  Atp  the  very  first  time 
Interval  selected.  Within  the  domain  of  these  limitations  we  have 
always  chosen  the  largest  time  steps  allowable. 

In  determining  the  time  step  for  a  new  time  line,  it  is  necessary 
to  scan  through  all  the  neighboring  da».a  points  on  the  initial  constant 
time  line  and  select  the  largest  time  interval  permissible.  A  data 
point  distribution,  along  a  constant  time  line,  which  makes  every  local 
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maximum  allowable  A  t  the  same  would  be  mop,t  desirable  for  the  sake  of 
computing  effort.  However,  generally  this  is  not  possible.  Neverthe¬ 
less,  with  this  in  mind,  we  occasionally  rearranged  the  distribution 
of  initial  data  points  in  the  course  of  the  computer  calculation. 


IV  RESULTS 


Due  to  the  complex  nature  of  the  flow  field,  the  calculation  is 
divided  into  three  stages,  as  shown  schematically  in  Fig,  6.  The 
first  stage  of  calculation  starts  from  t  **  0  and  ends  at  the  time 
when  the  head  of  the  rarefaction  reaches  the  center  of  the  sphere. 

The  second  stage  of  calculation  terminates  just  before  the  inward 
shock  reaches  the  center  of  the  sphere.  The  third  stage  includes  a 
triangular  region  bounded  by  the  constant  time  line  AB,  the  shock 
front  BD,  and  a  line  AD.  Point  A  is  placed  some  distance  away  from, 
the  center.  The  slope  of  AD  is  chosen  to  be  smaller  than  that  of  the 
local  I  characteristic,  (because  the  condition  a  t  <  thus  the 
reflected  second  shock  does  not  affect  the  region  to  the  right  of 
line  AD. 

In  order  to  make  the  results  more  general,  the  dependent  variables 
u,  c,  p  and  the  independent  variables  x  and  t  were  non-dimensionalized 
with  appropriate  parameters  of  the  atmospheric  air  and  the  initial 
total  energy  released  in  the  expanding  sphere.  The  non-dimensional 
variables  include: 

P 

pressure  -rr  =  — 

P, 

particle  velocity  a  = 

v  C, 


39  -  . 
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speed  of  sound  Q  =  — 

e 

density  ~  ~o~ 

» i 

radius  A  - 

1.  General  Features 

The  general  behavior  of  the  expansion  of  a  high  pressure  sphere 
can  best  be  shown  by  a  numerical  example.  Results  of  the  expansion  of 
a  high  pressure  sphere  initially  with  pressure  ratio  \  **  100,  density 
ratio  ,>|0*  1.16,  and  1.4  as  calculated  by  the  Hartree  scheme,  are 

shown  in  Fig.  7  to  Fig.  12. 

Fig.  7  shows  the  propagation  of  the  main  shock  followed  by  the 
contact  surface;  both  travel  outward  with  decreasing  speeds.  A  second 
shock  appears  at  a  later  time,  even  though  it  travels  to  the  left 
relative  to  the  gas,  for  a  moment  it  is  moving  with  an  absolute  velocity 
toward  the  right.  However,  as  the  second  shock  gains  strength  with 
time,  it  soon  acquires  an  inward  moving  speed  and  then  races  toward 
the  center  at  an  accelerated  pace. 

The  decay  of  the  peak  pressure,  7k  =  ,  where  p  is  the  pressure 

immediately  behind  the  main  shock,  is  shown  in  Fig.  8,  The  rate  of 
decay  decreases  as  the  main  shock  moving  to  a  larger  radius.  Pressure 
profiles  are  shown  in  Fig.  9  and  Fig.  10.  The  first  curve  (T*»  0)  in 
Fig.  9  is  the  initial  pressure  distribution,  and  the  second  curve 


-41- 


(t  =  0.0026)  shows  the  pressure  profile  shortly  after  the  release.  The 
rarefaction  wave  moves  toward  the  center  of  the  sphere  and  reaches  there 
at  t  -  0.023.  By  that  time  the  second  shock  is  already  formed.  Notice 
the  lowest  pressure  occurs  at  the  tail  cf  the  rarefaction  wave.  To  its 
right  pressure  increases  with  A,  the  non-dimensional  radius,  and  reaches 
a  peak  at  the  shock  front.  At  the  contact  surface,  which  is  marked  by 
a  bar,  pressure  is  continuous,  but  its  derivative  with  respect  to  A  is 
not.  For  clarity,  curves  shown  in  Fig.  10  are  drawn  with  a  different 
scale  from  that  used  in  Fig.  9.  If  this  is  unnoticed  the  general 
appearance  may  be  deceiving.  At  t  =  0.030  the  second  shock  already 
gained  considerable  strength,  and  attained  its  extreme  right  position. 

This  rapid  increase  in  shock  strength  is  primarily  due  to  the  diminishing 
of  pressure  in  front  of  the  second  shock.  After  the  main  shock  propagates 
quite  a  distance,  the  pressure  gradient  behind  it  gradually  decreases. 

This  can  be  seen  by  comparing  the  three  curves  on  the  extreme  right. 

Fig.  11  shows  the  profiles  of  the  particle  velocity  at  various  time 
instants.  Particles,  gaining  speed  through  the  expansion  wave,  reach 
their  peak  velocity  at  the  tail  of  the  rarefaction  wave.  It  is  interest¬ 
ing  to  see  that  at  x  =  0.045  the  particle  velocity  behind  the  second 
shock  becomes  negative,  which  means  gas  particles  are  pushed  back  toward 
the  center.  This  negative  velocity  reaches  very  high  value  just  before 
the  second  shock  implodes  on  the  center  of  the  sphere. 

c 

The  profiles  of  the  speed  of  sound  9  =  —  ,  are  presented  in  Fig. 

Cj 

12.  Since  for  a  polytropic  gas,  temperatures  is  proportional  to  square  of 
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the  velocity  of  sound,  therefore  these  curves  can  be  interpreted  as 
the  temperature  profile  with  square  root  of  temperature  ratio,  —  , 
as  ordinate.  Initially  the  high  pressure  gas  is  at  a  much  higher 
temperature  than  that  of  the  outside  air.  Immediately  after  the 
release,  gas  temperature  reduces  through  expansion  and  in  the  mean¬ 
time  air  temperature  rises  after  experiencing  a  shock.  However,  the 
difference  in  temperature  remains  quite  large  across  the  contact 
surface  which  separates  them.  It  is  shown  that  there  is  a  big  jump  in 
speed  of  sound  from  the  right  to  the  left  across  the  contact  surface. 

The  speed  of  sound  at  the  center  becomes  progressively  lower  as  a 
result  of  expansion,  but  the  implosion  of  the  second  shock  will  raise 
it  again. 

In  an  attempt  to  investigate  whether  the  distribution  of  the  flow 
properties  of  the  expanding  high  pressure  sphere  at  late  times  approaches 
that  of  Sedov's  blast  wave  solution,  the  distribution  of  the  three  flow 
properties  were  once  again  plotted  in  a  somewhat  different  fashion. 

Fig.  13  shows  the  ratio  of  local  pressure  to  the  corresponding  peak 
pressure  against  the  ratio  of  local  radius  to  the  corresponding  radius 
of  the  main  shock  front,  for  several  late  times.  The  solid  curve 
labeled  as  Sedov's  blast  wave  solution  was  based  on  the  numerical 
values  given  by  Sedov  (2).  Similar  curves  for  particle  velocity  and 
speed  of  sound  were  presented  in  Figures  14  and  15.  It  is  seen  that  the 
pressure  distribution  adjacent  to  the  main  shock  front  for  the  expanding 
sphere  approaches  that  of  the  Sedov's  blast  wave  solution  as  time 


-  43  - 


elapses.  However  this  conclusion  cannot  be  applied  to  the  other  two 
flow  properties,  particle  velocity  and  speed  of  sound.  In  fact  as  a 
whole  the  flow  field  resulted  from  the  expansion  of  a  high  pressure 
sphere  is  quite  different  from  that  of  Sedov's  blast  wave  solution. 
Nevertheless  the  destructive  effect  of  a  blast  wave  usually  is  deter¬ 
mined  by  the  peak  pressure  of  the  blast.  As  will  be  shown  later  in 
the  section  on  late-stage  equivalence,  the  peak  pressure  of  an 
expanding  high  pressure  sphere  can  be  approximated  by  the  corresponding 
Sedov's  solution  with  equal  total  initial  energy  in  a  region  when  the 
main  shock  has  propagated  some  distance,  and  the  peak  pressure  still 
above  10  atm.  Therefore,  for  a  rough  estimation  of  the  damaging  effect 
resulting  from  the  shock  pressure  of  the  expansion  of  a  high  pressure 
sphere,  the  corresponding  Sedov's  blast  wave  solution,  equal  energy,  can 
be  used. 

Some  numerical  results  of  the  expansion  of  high  pressure  spheres 
were  reported  previously  by  Brode  (26) .  Those  results  were  presented  in 
rather  sketchy  curves  with  ill-defined  boundaries.  Therefore,  it  is 
impossible  to  make  a  precise  assessment  of  their  accuracy.  However, 
judging  from  our  calculations  of  one  of  Brode' s  examples,  his  results 
are  in  considerable  error,  even  in  the  region  corresponding  to  our  first 
stage  of  calculation.  The  error  in  peak  pressure  is  estimated  to  be  of 
the  order  of  10%,  Even  larger  errors,  more  than  20%,  in  pressure 
occurred  in  the  rarefaction  wave  region.  Moreover,  the  pressure  profile 
3rode  identified  as  second  shock  was  so  widely  spread  that  it  could 
hardly  be  called  a  shock.  Nevertheless,  despite  all  fnese  discrepancies, 
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qualitatively  Brode's  results  are  in  general  agreement  with  our 
calculations. 

2.  Late-Stage  Equivalence 

A.  Late-stage  equivalence  in  hypervelocity  impact. 

Before  we  discuss  the  principle  of  late-stage.  equivalence  in 
expanding  high  pressure  spheres  it  seems  appropriate  to  review  some 
related  background  about  the  general  concept  of  late-stage  equivalence. 
The  principle  of  late-stage  equivalence  was  first  proposed  by  Walsh 
et  al  (27)  in  studying  effects  of  hypervelocity  impact.  It  stipulates 
that  projectiles  of  different  mass  and  velocity  striking  a  target  will 
give  rise  to  target  flows  very  nearly  identical  at  late  times,  pro¬ 
vided  the  product  of  mass  and  the  velocity  raised  to  the  d  power 
is  the  same  for  all  projectiles.  This  principle  is  assumed  to  be 
applicable  within  the  hydrodynamic  phase  of  the  hypervelocity  impact 
where  the  interaction  is  governed  by  equations  of  fluid  dynamics  prior 
to  the  onset  of  material  strength  effects.  If  two  impact  flows  during 
the  hydrodynamic  phase  of  the  interaction  are  equivalent,  then  the 
strength-dependent  phase  of  the  motions  will  also  be  equivalent.  There¬ 
fore,  the  property  of  late-stage  equivalence  enables  one  to  extrapolate 
the  existing  experimental  results  from  impacts  at  comparatively  low 
velocities  attained  in  controlled  experiments  to  the  meteoric  velocity 
which  is  not  accessible  to  laboratory  experiments. 

For  one -dimensional  impacts  of  materials  with  ideal  gas  equation 
of  state,  Walsh  et  al  (27)  showed  that  by  assuming  impacts  are  late- 
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stage  equivalent  on  the  basis  of  M0Ve  rather  than  dependent  on  Mq  and 
V0  individually  led  to  the  results  that  the  flow  is  of  self-similar 
type.  The  analytical  similarity  solution  and  a  solution  by  finite- 
difference  calculation  (SPUTTER  code)  of  an  ideal-gas  impact,  with 
?  "  1,4  were  found  to  be  in  good  agreement  at  late  times  on  the  basis 
of  equal  M0V0^‘^,  Recent  work  by  Chou  and  Bums  (28)  reported 
extensive  calculations  for  the  same  problem  employing  the  method  of 
characteristics  demonstrated  also  that  the  shock  fronts  produced  by 
different  one -dimensional  impacts  approach  each  other  in  position  and 
in  strength  at  late  times  provided  M0V0*,‘*  is  constant. 

In  the  case  of  axisymmetric  impact  of  ideal-gas  the  problem  is 
fully  characterized  by  given  values  of  Y ,  the  initial  density  f>t,  the 
characteristic  length  of  the  projectile  LQ,  and  the  impact  velocity  vQ, 
From  dimensional  considerations,  the  solution  must  be  functions  of  three 

T  7  t  V 

non-dimensional  variables,  —  .  and  V* ,  where  r  and  z  are 
cylindrical  coordinates  with  origin  at  the  center  of  the  initial  sur¬ 
face  of  impact.  Ine  assumption  of  late-stage  equivalence  on  the  basis 

-  4 

of  M,V,  or  l0V  is  to  require  that  the  solution  not  contain  LQ,  VQ 
as  individual  parameters.  The  solution  must  approach  at  late  time  a 
flow  expressible  in  the  form: 

/**'+*  Z^' 

?  ”  ?o  $1  ^  >  i*'  } 
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vhere  T  =  ,  R  ~  }  2  =  -j—  and  d' =  -j~.  This  represents  a 

similarity  solution;  i,e.,  the  number  of  independent  variables  is 
reduced  to  two  instead  of  the  original  three,  The  similarity  solution 

apparently  still  cannot  be  solved  analytically  without  additional 

y 

approximations.  However,  it  is  possible  from  the  form  of  the  solution 
to  determine  the  relation  which  must  exist  between  the  shock  decay  rate 

and  d’ t  shock  pressure  decays  as  2  *  >  and  between  the  rate  of 
increase  of  total  positive  axial  and  total  positive  radial  momenta  and 


d/.  These  integrated  momenta  increase  as  t 


.  This  is  because 


the  distances  increase  an  t  *  »  mass  engulfed  therefore  increases  as 

JsC  _ !_> 

t  >  and  velocities  decrease  as  X  ,  so  that  the  momentum 


within  corresponding  elements  of  volume  varies  as  t  .  These 
results  can  be  compared  with  .the  computed  solution  determined  without 
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assuming  late-stage  equivalence.  Walsh  (27)  calculated  a  problem  of 
right  circular  cylinder  of  unit  aspect  ratio  impacting  a  semi -infinite 
target  using  the  Eulerian  code,  and  found  that  the  calculated  solution 
agrees  at  late  times  with  the  similarity  solution.  The  value  of  ot' 
determined  from  the  shock  decay  rate  and  the  increase  of  momenta  is 

»  .59  which  corresponds  to  1.77.  These  results  are  for  Y  ■  1.5. 

For  solid-solid  impact,  Riney  (29)  proposed  a  vi3co-plastic  model 
which  bridges  the  transition  from  the  early  hydrodynamic  phase  to  the 
later  stages  when  strain-rate  and  strength  effects  become  important. 
Calculations  for  a  cylindrical  projectile  impacting  a  thick  target  of 
like  material  at  various  velocities  were  carried  out  to  the  point  where 
strain-rate  and  strength  effects  can  no  longer  be  neglected.  Upon  com¬ 
paring  the  flow  field  resulting  from  impacts  of  different  velocities 
Riney  found  that  the  flow  fields  are  essentially  equivalent  except  for 
the  initial  stage  of  the  cratering  process  if  the  geometrically  similar 
projectiles  are  chosen  to  have  the  same  kinetic  energy,  MQV0  ,  at  impact. 

In  short,  the  concept  of  late-stage  equivalence  of  hypervelocity 
intact  on  the  basis  of  M0\^*  is  established  on  numerical  evidence  rather 
than  on  theoretical  grounds.  However,  the  most  important  general  con¬ 
clusion  from  one-dimensional  impact,  the  axisymmetric  ideal  gas  impact 
and  solid-solid  impacts  is  that  impacts  are  equivalent  on  a  basis  inter¬ 
mediate  between  equal  projectile  momentum  and  equal  prdjectile  energy. 
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B.  Late-stage  equivalence  in  expanding  high  pressure  spheres 
In  view  of  the  evidence  of  the  late-stage  equivalence  in  hyper¬ 
velocity  impact,  it  leads  us  to  the  investigation  of  similar  equivalence 
in  expansion  of  high  pressure  spheres.  Based  on  the  numerical  results 
of  the  present  calculation,  late-stage  equivalence  is  found  to  exist  for 
the  expansion  of  high  pressure  spheres  with  various  initial  conditions, 

provided  the  initial  energy  released,  £c  =  3^ in  each  of  the 
% 

spheres  is  held  constant.  Late-stage  equivalence  is  assumed  to  exist 
if  the  peak  pressure  distribution  for  different  expanding  spheres  are 
the  same  at  late  times.  This  is  demonstrated  by  Fig.  16  where  peak 
pressure  vs.  radius  curves  of  expanding  spheres  with  equal  initial 
energy  EQ  but  different  initial  radius,  pressure  and  density  are  shown 
in  logarithmic  scales.  The  purpose  of  choosing  logarithmic  scales  is 
to  cover  the  extremely  wide  range  of  peak  pressures,  from  TT0  = 2000  to 
TTp  =  2  ,  in  a  single  plot.  Accompanying  these  curves  the  corresponding 
point  source  solution  is  also  included  in  the  plot.  The  point  source 
solution  represents  the  extreme  case  of  a  high  pressure  sphere  of  zero 
radius . 

In  the  high  pressure  region,  the  peak  pressure  distribution  of  the 
point  source  solution  can  be  derived  from  equations  given  by  Sedov  (2). 
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where  p  is  the  peak  pressure,  x  is  the  position  of  the  shock  front 

®  ft 

and  v-  ip  a  connt&nt:  proportional  to  E0,  tiie  tota?.  er-.^rgy  released 
initially.  According  to  the  point  t.crce  formulation  the  total  energy 
of  the  disturbed  gas  remains  constant.  Therefore,  ye  have 


X  % 

E0  »  J  *irx*Jx  *  l  7 $r 


The  first  term  is  the  kinetic  energy  and  the  second  term  is  the 
internal  energy.  Introducing  non-dimensional  quantities,  we  can  write 
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It  can  be  shown  that  the  sum  of  the  integrals  in  the  square  bracket  is 
a  function  of  JT  only  and  for  V  ■  1,4,  E  *  1.175E0.  Since  £  is  defined 


let 


be  the 


non-dimensional  radius  of  the  shock  front 


then  the  relation  between  peak  pressure  and  position  can  be  reduced  to 
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the  simple  form, 

1  ~3 

^  =  0.ISG1  As 

This  is  shown  by  the  straight  line  portion  of  the  modified  point  source 
solution  in  Fig.  13.  It  is  called  ''modified  point  source"  because 
when  the  peak  pressure  ratio  is  less_  than  100,  numerically  calculated 
results  are  continued  from  the  point  source  solution. 

For  peak  pressure  belcw  100,  the  point  source  solution  is 
modified  by  taking  into  consideration  the  ambient  pressure  p^.  This  is 
accomplished  by  carrying  out  a  computation  similar  to  the  calculation 
of  similarity  solution  of  blast  wave  by  standard  method  of  character¬ 
istics,  except  this  time  exact  shock  relations,  equations  (2-36)  to 
(2-38),  instead  of  the  strong  shock  approximations;,  equations  (2-39) 
to  (2-41) ,  are  utilized  in  the  calculation.  Numerical  results  obtained 
are  found  to  be  in  extremely  good  agreement  with  those  of  Goldstein's 
(3).  The  graphical  representation  of  the  modified  peak  pressure  dis¬ 
tribution  appears  in  Fig.  16  as  a  curved  line  with  a  positive  radius  of 
curvature . 

To  keep  the  initial  energy  released  constant,  spheres  with  higher 
initial  pressure  are  associated  with  smaller  radii.  It  is  seen  that 
the  peak  pressure  curve  of  a  sphere  with  higher  initial  pressure  decays 
faster  and  coalesces  with  those  started  at  lower  initial  pressure  and 
they  all  finally  approach  the  modified  point  source  solution.  Notice 
in  every  case,  immediately  after  the  release,  the  peak  pressure  is  much 
lower  than  the  corresponding  peak  pressure  of  the  point  source  solution, 
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but  it  decays  at  a  slower  rate  and  exceeds  the  modified  point  source 
solution  before  they  finally  coalesce. 

If  in  addition  we  require  the  initial  energy  released  to  be 

constant,  the  total  mass,  ^  poX*»  is  also  kept  constant,  the  late** 
stage  equivalence  exists  not  only  for  peak  pressure  of  the  main  shock, 
but  also  for  the  positions  for  both  the  main  shock  and  the  second  shock 
in  the  physical  plane  (i-X  plane).  This  is  revealed  in  Fig.  17  where 
the  main  shock  front,  contact  surface  and  the  second  shock  are  pre¬ 
sented  for  high  pressure  spheres  with  initial  pressure  ranging  from 
"IT,  *  50  to  TT0  «  500.  The  late-stage  equivalence  in  peak  pressure  can 
be  seen  in  Fig.  18,  This  time,  peak  pressure  curves  are  plotted  in 
linear  scales.  Perhaps  a  presentation  like  this  can  be  better 
appreciated.  Once  again,  the  four  curves,  for  TTt  «  500,  % «  200, 

TT,  «•  100  and  IT,  ■  50  converge  together  and  approach  the  modified  point 
source  solution  at  late  times. 
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V  CONCLUSION 


The  Hartree  scheme  was  successfully  applied  to  the  solution  of 
the  flow  field  of  the  expansion  of  high  pressure  spheres.  A  special 
technique  designed  to  handle  the  starting  singularity  was  found  to  be 
accurate.  Results  obtained  by  this  method  for  the  initial  behavior 
of  the  spherical  expansion  agree  very  well  with  those  calculated  by 
McFadden's  analytical  solution. 

Tiie  formation  of  a  second  shock  from  the  previously  continuous 
flow  field  was  found  to  exist  at  the  tail  of  the  left  traveling  rare¬ 
faction  wave.  The  strength  of  the  second  shock  remains  rather  small 
before  the  head  of  the  rarefaction  wave  reaches  the  center  of  the 
sphere,  then  it  grows  rapidly  to  very  high  strength.  The  path  of  the 
second  shock  in  the  physical  plane  is  accurately  defined  and  some 
interesting  features  are  revealed. 

The  numerical  results  of  a  typical  example  of  the  expansion  of  a 
high  pressure  sphere  were  presented  in  graphs  showing  the  propagation 
of  the  main  6hock,  the  development  of  the  second  shock  and  the  propaga¬ 
tion  of  the  contact  surface  which  separates  the  two  fluids.  The 
boundaries  of  these  discontinuity  surfaces  were  accurately  revealed. 
Profiles  of  pressure,  particle  velocity  and  speed  of  sound  at  different 
times  after  the  release  were  also  depicted. 

Based  on  numerous  calculations,  it  is  shown  that  "late-stage 
equivalence"  exists  in  the  expansion  of  high  pressure  spheres  into  the 
atmosphere,  provided  the  initial  total  energy  in  each  of  the  spheres  is 
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kept  constant.  If  the  initial  condition  is  further  limited  to  equal 
masses,  then  the  late -stage  equivalence  exists  not  only  for  peak 
pressure  of  the  main  shock  but  also  for  the  position  of  both  the  min 
and  the  second  shock  in  the  x-t  plane. 

The  method  developed  in  this  report  can  be  extended  to  solve 
problems  for  perfect  fluids  with  equations  of  state  other  than  that 
of  the  poly tropic  gas.  It  may  also  be  extended  to  radially  symmetric 
problems  with  non-uniform  initial  properties  nnd  non-zero  initial 
particle  velocities. 
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Figure  4  Strength  of  Second  Shock  vs  Time  for 

Different  Initiation  Conditions  (  100,  V*  1.16,  1.4) 


Schematic  of  Reg: 


Showing  the  Main  Shock,  Contact  Surface  and  the 


Pressure  (  vs  Radius  C  A,=  f 100,  7,.  1.16,  1.4) 


Figure  9  Pressure  (  W-  ■£  )  vs  Radius  (A=-<r)  at  Indicated  Times  (  7J »  100,  1.16,  1.4) 
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Figure  16  Peak  Pressure  (  JJjs|[)  vs  Radius  (  Aj2^)  for 

Expanding  Spheres  with  Equal  Initial  Energy,  but 
Different  Initial  Pressure  Density  and  Radius 


ressure  (J7j=y)  vs  Radius  C^j  =  ^)  for  Expanding  Spheres  with 
Initial  Energy  and  Equal  Initial  Mass 
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